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Abstract 



The presented article contains a 3D mesh generation routine optimized with the 
Metropohs algorithm. The procedure enables to produce meshes of a prescribed volume 
Vq of elements. The finite volume meshes are used with the Finite Element approach. 
The FEM analysis enables to deal with a set of coupled nonlinear differential equations 
that describes the electrodiffusional problem. Mesh quality and accuracy of FEM solu- 
tions are also examined. High quality of FEM type space-dependent approximation and 
correctness of discrete approximation in time are ensured by finding solutions to the 3D 
Laplace problem and to the 3D diffusion equation, respectively. Their comparison with 
analytical solutions confirms accuracy of obtained approximations. 



1 Introduction 



One from the most important physical processes is electrodiffusion. It describes both 
diffusional motion of mass and charge flow due to applied electric field. The electric po- 
tential distribution is govern by the Poisson equation and total transport of particles is 
given in terms of the continuity equation [1] . The significance of this equation is broadly 
described in existing physical, chemical and biological literature p] and lots of scientific 
articles, particularly, those which concern properties of nano and micro transport. 
Mathematically, equations of electrodiffusion constitute a set of coupled nonlinear equa- 
tions where the Laplace operator [21 SI E] appears together with the first order partial 
time derivative. The Laplace operator is the basic operator met in many physical situa- 
tions. Thus the first step to deal with the electrodiffusional problem is to approximate 
the solution of the Laplace equation with help of the Finite Elements Method. Practi- 
cally, it means that an appropriate mesh should be designed for a prescribed 3D domain. 
The mesh must fit well to the physical conditions like e. g. symmetry of the problem. 
Therefore, different mesh shapes could be desired (spherical, cylindrical, conical or cubic) 
up to problem. After having accurate basic spatial solutions on appropriate meshes, the 
problem should be extended to the time-dependent case of the diffusion equation by find- 
ing discrete approximation in time. It could be done by means of truncated Taylor series 
or other single-step procedures like the Crank-Nicolson scheme [HI E] or the Gurtin's 
approach to finite element approximation in terms of variational principle [8]. From now, 
further extension of above-presented computations involving non-linear terms could be 
easily implemented and numerically solved using the Newton's method |9]. 



2 Equation of electrodiffusion 

The equation of electrodiffusion [1] has the form 



Ot 



where rij is a number of z-th ions, - electric potential, Di - the diffusion coefficient of 
i-th particles, ks - Boltzmann constant, T - temperature, Zi - valence of the i-th kind of 
ions, e - electric charge. To find the electric potential the Poisson equation must be 
solved 

V'^(eoeV0) + p = O (2) 

where p = '^^Pi, pi = Zicrii. Thus equations ([I]) and (El) both constitute the system of 
coupled equations. 

2.1 FEM approach 

Next, they can be solved numerically using the Finite Element Method [6] where the 
problem is represented as 

[ v^A(u)dn = / [viAi{u) + V2A2{u) + V3A3{u)]dn = 



j v^B(u)dr = j [5i5i(u) + v^B^in) + v^B^{n)]dT = 



(3) 



where 
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n_ 




(4) 



and V and v are sets of arbitrary functions equal in number to the number of equations (or 
components of u) involved. Ai{u),A2{u) and ^3(u) are given by the following formulas, 
respectively 
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+ [0, 1, 0]-u (5) 



^3(u) = eoeV'^Vu^ 
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where ki — , i — {+, —}. An expression B(u) gives the boundary conditions on F, 

ksT 

however, we choose a forced type of boundary conditions on F i. e. 

- 06o«n = 
l^i ^i, boun 0. 



(6) 



Let us substitute u ~ A^^^I 
In that way, we end up with the Galer' 



= Nu and put v = ^^w^Kub where Wb = N^. 
dn formulation of the problem. 



2.1.1 Discrete approximation in time 

In turn, we can approximate the nodal electric potential 

(j){x, y, z,t) = ^ Na{x, y, z)<p''{t) 



and number of particles 



at a time tn by 



n, 



{x, y, z,t) = ^ Na{x, y, z)n'}{t) 



and taking advantage from an evaluation of u„ = 



0n 
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(7) 
(8) 

(9) 
(10) 

in the Taylor series we obtain 

(11) 



where /3 takes values from [0, 1] and At denotes time step. After incorporating it into a 
general form of time-dependent equations ^{1,2} (u) 

K(u) + Cii = (12) 

where K{u) represents these parts of /l{i^2}(u) with a space-dependent operator we get 
time approximation for a given node 

When /3 = 1 an approximate solution to the semi-discrete equations at each time t„ is 
given by the Euler „backward" scheme 

K{Nu) + CN^u„ = CN^u„_i (14) 

otherwise, the expression for is as follows 

^^^^Wt'''' ^ ~ ^^"^^"""-i + ^K{NJnU) = 0. (15) 

2.1.2 Space— dependent term 

After integration A'^ii'(Nu) by parts, one obtains a mixed set of hnear and nonlinear 
equations 

+it L ^ ^« " ^^-^^ ^ ° 

+Xtl^^'Y. ^'^ (^n'" - n-L\) dn = 

-J {WNbfeoeV (^N^r^ dft + z^e J N' ^X^iVan+'"j dfi 



dn^o 



6 = 1,...,M 

(16) 



3 



and a corresponding set of boundary terms for and where i = {+, — } 
+eeo^iV^|^|v ||^iV,rj|rfr = 6=1,. ..,M 

(17) 

d 

where — denotes derivative normal to T. Presented-above spatially temporal discretiza- 
on 

tion has been done for the case with /3 = 1 at each node. If the forced boundary conditions 
(see Eq. (Q) are imposed on and r„, respectively, then all terms in Eq. (fT7|) can be 
neglected by restricting the choice of A*";, functions to those which equal on F. Let's 
denote integrals from Eq. ( !TB|) as 

Kac= [ {VNtfNaVNjn = J2 [ {VN,fNaVNjn'' = J2Kc 

= I iSN.Y VNjn = J2[ (ViV,)^ vKdn^ = K'f (is) 

Jn g Jn- 

Kl= f NaN^dn = V / NaN^dn" = V K^f. 
Jn g Jn^ 

where ^g with e = 1, . . . , E denotes sum over elements. In 3D space tetrahedral elements 
seem to be a natural choice of finite volume elements. Then indices a, b, c take four values 
each (an element has four nodes) from the set of 1, . . . , M values. 

2.2 Tetrahedral elements 

For tetrahedral linear elements shape functions Ni can be assumed as equal area coordi- 
nates Li given by the formula 

U = , « = 1,2,3,4 (19) 

where is a volume of tetrahedron. The following integration formula can be useful 



I ^L^.L^LUxdydz = ^lEll^ Qye 

Jn^ ' ' ' ^ (« + /3 + /i + z/ + 3)! 



(20) 
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in order to calculate integrals K^^, K^'^ and K^'^. Shape functions for linear elements 
are Na = La for a = 1, 2, 3, 4. This gives 

^ {h%^ + C^Ce + Sd^ 



mv 



k'f = [ {VL'fVLjn' = {h%a + c'Ca + d'da) (21) 

K^f = / L^L^dl]'^ = — - when a 7^ b 
Jne 5! 
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Finally, we have 

^ V e / e e ' 

W ( E ) + ^- E ^''^ + E ^'4^"'" = 

^ \ e / e e 

eoer E = l^^l^ E ^ (^^'"^ - ^"'1 

e e 

where a,b,c = 1, . . . ,M 



(22) 



where /*'^ = — ^^^-^^^'"^^ ^n-i with z = +, — and only these nodes a,b and c that 

participate in the particular element e can give a non-zero contribution to the sums of 
the general type ^g-ft'^- 

Let's assume that all values of are known at a time t„. Then we can solve the third 
equation in ([22D obtaining the result 
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{E^'"} (^E^'") (^^"-^"") (23) 



where j^^^-^a^j denotes elements of the inverse matrix to K. After substitution the 
solution for to Eq. (!22|) we get 
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where z = +, — . 



2.2.1 Newton's method 

The above written set of equations is of a nonhnear type. Let's denote all of them as 

(25) 



Thus to solve it we have to employ the iterative Newton's method J9j. It means that 
we have to start from an initial guess of {tzq}^^ values. And during next iterations for 



k = 0,l,... 



(26) 



the solution should be achieved. F' (n\, . . . ,n^^) denotes the following matrix of partial 
derivatives: 



dF 



M 



dF^ 



dF^ 



where 



dF'' 
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3 Diffusion equation 

Putting A; = in the equation of electrodiffusion we neglect the electrostatic term. It 
leads to the following equation describing diffusion in JJ" [H [3] 



Ut-DAu = in 3fJ« x (0,oo), 
u = g on 3?" X {t = 0}, 



(29) 



where D denotes a diffusion coefficient. This kind of equation represents an initial value 
problem. Assuming that the considered domain Q in 3?^ is of a cubic type [xmin,Xmax] x 
[ymin,ymax] X [zmin, Zmax] let US take M = as a bouudary condition. Now we seek 
a solution of the equation ( l29i) which satisfies this boundary condition and prescribed 
initial condition at the time t = 0. The solution of the equation is approximated by the 
triple sum 



oo oo oo 



^) = E E E ^o,fe.,fc„fc.e- 



{kl+kl+kl)Dt 



sin{kxx) sin{kyy) sm.{kzz), (30) 



where VQ^kx^ky,k^ are unknown coefficients that must be determinated from the initial con- 
dition: 



oo oo oo 



5, = m(-, 0) = E E E ^O'l^-^'^yk. sin(/c^x) sm{kyy) sin(A;^z). 



(31) 
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In the case of the domain being [0 tt] x [0 tt] x [0 tt] and g = const the solution has the 
form 

64 °° 1 

u(x,y,z,t) = ^ e~^^'+*^^+^'^^*sin(A:^3:)sin(A;„y)sin(fc^^). (32) 

kx t^z — 1 1 3 5 . . . 

In the case of cyhndrical domain defined by r e [0,ro], 9 G [0,27r) and z G [0,7r], and 
with the boundary condition of the form u = the solution of Eq. fl29|) can be expanded 
in an absolutely and uniformly convergent series of the form 

oo oo 

M(r, e,z,t) = Y^ JniK,mr/ro) (a„,™ cos(n^) + bn,m sin(n^)) sin(A;,z)e-«'="-'"/"°)'+'^')^* 

n=0 m,kz=l 

(33) 

where fc„^m are the zeros of the Bessel functions and a„^m, bn,m are constants that must 
be found from the initial condition by making use of the orthogonality relation for the 
trigonometric and Bessel functions. 



4 Laplace equation 

On the other hand, assuming that the time derivative in the diffusion equation ( l29l) equals 
and putting D = 1 we end up with the boundary value problem of the Laplace type 

[HE] 

Am = in fi, ^^^^ 



u = g on dfl. 

Let's consider Q being a cubic domain i. e. [0 tt] x [0 tt] x [0 tt]. And for the g function 
equals everything on dQ apart from g{x = n, y, z) = (po one can approximate the exact 
solution by 



1600 sinh (y/n"^ + m^x) sin (ny) sin (mz) 
(x,y,z) = — — > , (35) 



vr^ nm sinh (\/n^ + m^vr) 

n,m,=l,3,.-- ^ ' 

For (j){x,y,z) = v{r) where r = (x^ + + z'^Y^'^ the Laplace equation has the solution 
defined in 3?^ for r ^ 

'^^^'^'^^ = 3«(3) {x^ + y'+z^y/^ ^''^ 
where a{3) denotes volume of B(0,1) in and equals 



r(3/2 + l)' 

5 Three-dimensional mesh generation 

Below are listed a few technical remarks referring to the mesh generation routine applied 
to obtain a designed 3D mesh. 
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5.1 Initial mesh 



An initial mesh is built on the basis of main surface nodes {outer nodes) which define 
a figure's shape. The whole figure is considered as divided into perpendicular to ^-axis 
layers. Thus the outer nodes are distributed on the edges of layers. In the center of 
each layer and also in the middle between two layers are located inner nodes. They are 
connected with outer nodes creating in this way the main figure's construction. Initial 
mesh elements obtained in such a manner are of tetrahedral shape. 

5.2 Figure's surface 

The boundary of the figure is defined by set of surface equations for vertical and horizontal 
segment lines linking outer nodes. After each mesh iteration new nodes are created and 
labeled as outer or inner ones according to surface equations. Moreover, the location of 
each node (i. e. on which exactly vertical, horizontal line or surface patch the node is 
lying) is also stored. 

5.3 New elements creation 

New elements are created by a division of already existing elements. At the beginning 
of the routine, the surface of division mainly connects a new node born on the longest 
element edge with two other nodes belonging to that mesh element and one node from 
the divided edge. The procedure constitutes a 3D extension of the 2D mesh generation 
routine described already in [H]. However, during the routine a number of small elements 
is increasing, and the division of the longest bar is not anymore the optimal way of 
proceeding. That is why, before choosing an edge to the division volume of elements 
common to it is checked. The edge that will not produce new elements having its volume 
smaller than an assumed critical volume is chosen to be cut. 

5.4 Mesh optimization — Metropolis algorithm 

The optimization is done with help of the Metropolis algorithm. The system energy is 
calculated as a sum of discrepancies between an element volume and assumed element 
volume Vq = where Hq denotes a prescribed length of the edge 



Thus the smaller a degeneracy from a designed volume distribution the more optimal 
state. The Metropolis routine starts from a nodal configuration given by described above 
procedure. The main point is to reach the optimal global configuration by ascertaining 
local optimal states. They arise from such a configuration of i-th node and its neighboring 
nodes which gives smaller energy E'^. This partial energy is calculated from the sum 
Eq. fl371) taken over elements containing the node of interest. To compute new positions 
for each node (giving new configuration) the following expression is put forward 




(37) 



e 






(39) 
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Figure 1: The figure presents four regular sliapes (cube, cylinder, sphere and cone) ob- 
tained with the non-optimized routine. Colorbars show variations in final elements vol- 
ume. 
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Figure 2: The figure presents a distribution of a normalized mesli edge lengtli L/ho 
for meshes in a cubic domain [0,7r] x [0, tt] x [0,7r] with a) unique elements created for 
ho = 0.22 and b) meshes optimized with the Metropolis algorithm with ho = 0.2. 



where kg denotes a shifting strength and Pi — Pj is the length of ij edge. The value 
of kg determines the strength of a nodal shift and varies from to 1. It is also worthy 
considering to choose its value as a random number from uniform distribution U{0, 1). 
Within the Metropolis routine the transition probability is calculated by the formula 

P{P^ ^ P^,nen.) = e-^""^/'-^ (40) 

where ks is a Boltzmann constant (here set as 1), T temperature and AEi = Ei^new — Ei 
is a difference between energies of these two states. If a value of P is greater than a 
random number from U{0, 1) new state is accepted. Otherwise the old one is preserved. 
All above-described local Metropolis steps can lead to different global configurations. 
Therefore, for each division number this distribution of nodes which gives a lower energy 
of the total system should be kept. To find this, again the Metropolis rule is employed, 
but this time changes in the total energy of the whole system are examined. To estimate 
maximal temperature T^ax (in expression ( HOl) ) a range of changes in potential energy 
corresponding to current number of elements must be found. Moreover, in each global 
Metropolis step temperature might decrease according to T = riT where parameter rj < 1. 



5.5 Delaunay reconfiguration routine 

To improve mesh quality the following transformations are applied |6]. 

• Three elements common to an edge are transformed to two elements, when one of 
the elements fails to satisfy the Delaunay criterion [TCTj . 

• Four elements common to an edge are transformed to a new configuration of four 
elements, when one of the elements does not meet the Delaunay criterion [10]. The 
new pattern is chosen from two different possibilities [H]. 

Additionally, too small boundary elements could be destructed by a projection of its 
internal node to the center of the outer patch of element that is opposite to it. Such 
approach is justified in the case of boundary elements, however, in the case of internal 
ones leads to creation of so-called irregular nodes. 
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Figure 3: The picture presents a distribution of a normalized mesh element volume V/Vq 
for the same domain as in Fig.[2]and for a) unique elements; b) elements with Vq = 9 x 10^^ 
optimized with the Metropolis algorithm. 

6 Results 



6.1 The Laplace equation 

The accuracy of FEM approximation of the Laplace equation on different meshes were 
examined. Numerical results vs. analytical ones for cubic and spherical domain are 
presented in Fig. (jlj). The relative difference between both analytical and numerical 
solutions has been calculated as 



(panal)/ max(0, 



anal i 



(41) 



The Laplace equation has been solved for the cubic domain [0, vr] x [0, vr] x [0, vr] with 
potential function = everywhere on the boundary F apart from one its side at x = vr 
where potential = 1 and for the spherical domain with the boundary conditions imposed 
by putting an elementary charge outside the sphere in [0,0,27r]. The exact solutions for 
both considered cases are evaluated precisely in Sec. 4. FEM approximation has been 
computed for the ,, linear" order of tetrahedron [6]. However, the comparison between 
both orders of approximation i .e ,, linear" and ,, quadratic" for the uniform mesh with 
Kzem = 0.0202 has been performed. The formulas for higher orders of approximation i. 
e. quadratic and cubic can be found in [B]. The results show that mean discrepancy be- 
tween numerical and analytical solutions calculated according to Eq. fHT]) for the Laplace 
equation equal —0.0061 ± 0.0153 in the linear case and 0.0004 ± 0.0082 in quadratic ap- 
proximation, respectively. Thus in further studies linear approximation will be used as 
being sufficiently accurate. In the case of cubic domain a mesh of unique element volume 
(non-optimized) has been applied in contrast to the spherical domain where mesh has 
been used after its enhancement with both the Metropolis algorithm and the Delaunay 
routine. 



6.2 The diffusion equation 

To test accuracy of discrete approximation in time the equation of diffusion (see Sec. 3) 
has been examined. To find P{x, y, z, t) particles distribution at tn moment in time the 
Taylor expansion (see Eq. [TTj) with /3 = has been applied. The diffusion equation with 
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Figure 4: The picture shows a) FEM approximation (hnear) of the Laplace problem at the 
center of the cubic domain i. e. at z = 7r/2 (as described in Sec. 4) vs. the analytical result 
together with b) a distribution of differences between numerical and analytical solutions 
obtained for each node in Q domain; c) FEM approximation (linear) vs. the exact solution 
of the Laplace equation for the spherical domain with the boundary values defined by 
putting an elementary charge outside the sphere in [0, 0, 27i] b) the volume distribution 
in the spherical domain optimized with the Metropolis algorithm with elements of a 
prescribed volume Vq equals 0.0075. 
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Figure 5: The picture shows a) FEM approximation (hnear) of the initial-boundary value 
problem in the center of the cubic domain i. e. at z = 7r/2 and at the time t = 0.19 with 
At = 0.01 [units] vs. the analytical result together with b) a distribution of differences 
(Eq. |4T|) between numerical and analytical solutions obtained for each node in Q domain; 
c) FEM approximation (linear) of the diffusion equation for the cylindrical domain at 
z = 7r/2 and at the following times: to = (blue), tmid = 0.5 (black) and t^nd = 1-0 
(red) with At = 0.01 [units]; d) volume profile of elements within the cylindrical domain 
where Vq = 0.015 and ho = 0.5; mesh quality were enhanced with help of the Metropolis 
algorithm. 



the following initial condition g{-, 0) = —x{x — T^)y{y — tt)z{z — vr) for the cubic domain 
and g{-, 0) = |(r — Ro)z{z — tt)\ where Rq denotes a radius of the cylindrical domain has 
been solved. The boundary value of the P{x, y, z, t) is set as 0. Results for both domains: 
cubic (with uniform elements - see Fig. [3K) and cylindrical (see Fig. [1]), this one tuned 
to the designed element volume with the Metropolis recipe, are shown in Fig. [5l In the 
case of cylindrical domain mean values of the ratio P(x, ti)/P{x, tj+io) calculated at each 
point of domain for times i = 0,10, ... , 100 have the average value 1.2699 ± 0.0049. 



6.3 The electrodiffusion equation 

The system of coupled equations describing process of electrodiffusion ([T]) written in terms 
of the FEM method (see Eq. ([22])) with the following values of constants: = = 
0.05, = D_ = 0.05 and the time step equals 0.01 has been numerically solved by 
using the Newton's algorithm. The boundary values of ?T,+ ,n_,0 are set as 1 at a; = 0, 
X = n, y = 0, y = n, and z = and equal 2 at z = n. An initial guess of n^,n^ and (p 
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Figure 6: The picture shows FEM approximation (hnear) of the system of electrodiffusion 
equations being the boundary value problem. The solutions are depicted in the center of 
the cubic domain i. e. ai z = -k /2 and at the time T = 0.39 with the values of parameters 
At = 0.01 [units], = k_ = 0.05, = D_ = 0.05 [units] a) the distribution of cations 
n+; b) the profile of the potential 0; computations have been performed on the uniform 
mesh with the volume of tetrahedron Vq = 0.01 and the element size ho = 0.44. 

distributions has been chosen as everywhere in the domain apart from its boundaries. 
The system of equations has been computed up to the final time T = 0.39. Fig. [6]presents 
obtained profiles of cations n+ and the potential (p at the center of the domain i. e. at 
z = tt/2. There is no visible difference between cations n_|_ and anions n_ distributions 
so the latter is not shown. The maximum of — n+^j is decaying from 0.023 for i = 

[t = 0) to 0.0093 for i = 37 {t = 0.37). The maximum of difference between and 
n_ distributions computed at each node equals 1.3e — 09. Employing physical notion it 
means that the system of charged particles is electroneutral. Moreover, the system of 
particles is tending to its stationary state. 

Additionally, components jx,jy and jz of the total flux of n+ particles flowing through 
the domain Q have been computed. They are shown in Fig. [71 Presence of a difference 
in an amount of n+ particles at the both sides of Q in the 2;-direction i. e. at z = and 
z = n causes a non-zero flow along z axis whereas a lack of such a difference in the two 
other directions i. e. x and y leads to the vanishing fiowsja; and jy in the center of the 
domain. 



7 Conclusions 

The presented software offers a 3D mesh generation routine as well as its further appli- 
cation to the 3D electrodiffusional problem. 

The proposed mesh generator offers a confident way to creature a quite uniform mesh 
built with elements having desired volume. Mesh elements have been adjusted to as- 
sumed sizes by making use of both the Metropolis algorithm and the Delaunay criterion. 
Mesh quality depicted in histograms occurs to be fairly satisfactory. Moreover, goodness 
of obtained meshes together with robustness of their applications to the Finite Element 
Method have been also tested by solving the 3D Laplace problem and the 3D diffusion 
equation on them. Comparison between these numerical solutions and analytical results 
shows very good agreement. 

To find solutions to a nonlinear problem defined by a system of coupled equations de- 



14 




Figure 7: The picture presents FEM approximation (linear) to the flux components ob- 
tained for the electrodiffusional problem. The solutions are shown at the center of the 
cubic domain i. e. at z — 7r/2 and at the time T = 0.39. Computations have been done 
with the following values of parameters: the time step At = 0.01 [units], nonlinear term 
multipliers kj^ = = 0.05 [units], diffusion coefficients = = 0.05 [units]. This 
figure shows approximated solutions for flux components for cations n+ in the y direction 
a) and in the z direction b) . The computations have been performed on the uniform mesh 
with the volume of tetrahedron Vq — 0.01 and the element size ho — 0.44. 

scribing electrodiffusion the FEM approach and the Newton method have been jointly 
applied. Analysis of obtained results conflrms usefulness of the presented solver to deal 
with nonlinear differential problems. 
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